Molecular dynamics simulations of the contact angle between water droplets and graphite surfaces 
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Wetting is a widespread phenomenon, most prominent in a number of cases, both in nature and technology. 
Droplets of pure water with initial radius ranging from 20 to 80 [A] spreading on graphitic surfaces are studied 
by molecular dynamics simulations. The equilibrium contact angle is determined and the transition to the 
macroscopic limit is discussed using Young equation in its modified form. While the largest droplets are almost 
perfectly spherical, the profiles of the smallest ones are no more properly described by a circle. For the sake of 
accuracy, we employ a more general fitting procedure based on local averages. Furthermore, our results reveal 
that there is a possible transition to the macroscopic limit. The modified Young equation is particularly precise 
for characteristic lengths (radii and contact-line curvatures) around 40 [A]. 
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I. INTRODUCTION 

Wettability is a long-standing issue primarily addressed by 
looking at the angle of contact at the edge of the interface 
between a liquid and a solid (sessile droplet method) QJ. In 
spite of the simplicity of the formulation of the problem, this 
procedure is at the basis of many investigations devised to as- 
sess the behavior of a liquid or a material in a number of in- 
dustrial processes and applications (see Ref. [2] for a review), 
and it has been demanding valuable efforts, both experimental 
and theoretical [3-6]. Computer simulations provide useful 
guidelines for their power to deal with the complexity of large 
assemblies of interacting components. In classical molecular 
dynamics studies, the systems are described at the molecular 
level according to the laws of classical mechanics and electro- 
dynamics. For these reasons, this approach proves to be both 
versatile and accurate. Recent breakthroughs have triggered a 
burst of interest in the properties of graphene GH9]. Promising 
applications in a variety of fields could indeed develop. Our 
interest in graphitic materials is related to their optimal disper- 
sion in polymeric matrices for the best manufacturing of com- 
posite materials lfT0WT2l . In that respect, the wetting prop- 
erties of graphitic surfaces by water are of course the starting 
point of any subsequent investigation. In particular, the transi- 
tion to the macroscopic limit is essential for any further study 
adopting more advanced modelization schemes |[T3l[T4l . Our 
analysis of profiles differs from the well-established method 
|[T5l[T6ll in that we approximate them by a more general curve 
than a circle. This way of proceeding is especially neces- 
sary for small droplets, exhibiting the largest deviations from 
the predictions of Young equation and from a spherical cap. 
This last aspect has already been recognized in previous stud- 
ies ifTTll . Yet, we also propose an analytic expression for the 
oxygen-oxygen radial distribution function of water. 



II. SIMULATIONS 

All molecular dynamics simulations are performed with 
LAMMPS 1 18 1, a code that supports parallelization optimally 
lfT9l . Numerical integration is accomplished with the algo- 
rithm rRESPA, allowing to deal with multiple time step sizes 



lfT8l . We choose a time step of 2 [fs] for non-bonded interac- 
tions and of 1 [fs] for bonded interactions. All Lennard- Jones 
forces are evaluated with a potential of the type 12-6. The in- 
built CHARMM force field l20l is employed to prevent van 
der Waals interactions from decaying abruptly at the cutoff 
distance: the forces are smoothly corrected to zero from 10 
to 12 [A]. Pairwise Coulomb interactions within a distance of 
12 [A] are calculated in the real space and beyond this value 
with a particle-particle, particle-mesh method; the precision is 
set to 10~ 4 . The pairwise interactions among atoms separated 
by one or more bonds are neglected. The neighbor lists are 
always updated at every time step. These general settings are 
always applied, unless specified otherwise. 



A. Water 

Throughout our work we use for the water the SPC/Fw 
model introduced in Ref. ETI . We address the reader to 
this detailed study for the definitions (partial charges, equi- 
librium distances, interaction parameters, etc.). We start from 
512 molecules of water arranged regularly in a cubic box of 
side length 30 [A] with periodic boundary conditions. We let 
the system evolve for 200 [ps] in the ensemble NPT (Nose- 
Hoover integration). This simulation is performed as equili- 
bration with a single time step size of 1 [fs]. The target tem- 
perature and pressure are 298.16 [K] and 1 [atm], respectively. 
Since the main purpose is to reproduce experimental densities, 
we choose the parameters that control the convergence so as 
to fix the temperature and the volume, while we still let the 
pressure fluctuate. It is in fact well-known that the pressure is 
a very sensitive function of the volume and difficult to equili- 
brate accurately [22]. We then let the system evolve for other 
0.5 [ns] at NVE conditions; the final configuration of this dy- 
namics is replicated and used to obtain the droplets of water. 



B. Droplets and graphitic substrate 

All starting configurations are formed by two parallel 
planes of graphene and a semisphere of molecules of wa- 
ter. The planes of graphene are separated by 3.4 [A] with 
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the lower one translated by the vector (l/2,y3l/2) with re- 
spect to the first; / = 1.42 [A] is the length of the bonds 
among the carbon atoms. The two planes are approximately 
squared with side of at least 30 [A] larger than the diam- 
eter of the semisphere. The semisphere is centered above 
the planes of graphene at a distance of 3 [A] from the up- 
per plane. The boundaries are periodic and two images of the 
droplets are separated by at least 100 [A] in the z direction. 
The mass of the carbon atoms is mc = 12.011 [g/mol]. The 
water and the graphene interact via van der Waals forces be- 
tween the atoms of carbon and oxygen with force field param- 
eters £ CO = 0.0478 [kcal/mol] and ct C o = 3.581 [A] E3ll24l . 
The system is equilibrated for 0.5 [ns] in the ensemble NVT 
(Nose-Hoover thermostat) with the temperature of the water 
maintained at 298.16 [K]. The system is studied during a fur- 
ther evolution of 1 [ns] in the microcanonical ensemble by 
gathering frames at every 0.5 [ps]. 



III. ANALYSIS 
A. Water 

The density of water is calculated using the formula p = 
N(m + 2m H )/(3 • 0.602 • Vdom)- N is the total number of 
atoms and Vd om is the volume in A 3 of the cubic domain re- 
sulting from the preliminary simulation of equilibration. By 
using this formula the result is expressed in g/cm 3 . The 
oxygen-oxygen radial distribution function g(r) is defined by 
(N/3V dom )g(r)4nr 2 Ar = W(r). W(r) is the average number 
of oxygen atoms in a shell of width Ar at a distance r from a 
given oxygen atom. We choose Ar = 0.05 [A] and of course 
g{r) is unitless. We also calculate the cumulative probability 
P(r) = (3/N) Jo W(r)nAr, which yields the fraction of oxy- 
gen atoms within a distance r from a given atom of oxygen. 
n = 1/Ar is the number of shells per unit length. (By writing 
the above integral as a discrete sum, that sum would run over 
the number of shells.) 



B. Droplets 

The contact angle 9 is defined by the tangent at the con- 
tact line, the edge of the interface between the solid and liq- 
uid phases (see Fig. B). The contact angle of a macroscopic 
droplet with spherical symmetry is well described by the equa- 
tion cos0 = (y sv — 7si)/Tlv- The y's are the surface/interfacial 
tensions. The subscripts s, 1 and v stand for solid, liquid 
and vapor, respectively. The above relation is generally re- 
ferred to as Young equation [ 1 ]. Especially for small droplets, 
Young equation needs a corrective term that accounts for the 
line tension, leading to cos0 = cos(L — kKj^R). As the 
notation suggests, cos 9„ = (y sv — 7si) / Hv is Young equation, 
which yields the contact angle &*> in the macroscopic limit. 
R is the radius (or curvature) of the contact line and K is 
the line tension. If 9 < 90°, we speak about hydrophilic 
behavior; hydrophobic if 9 > 90°. Complete wetting cor- 
responds to 0=0°. The discussion of the contact angle is 



^dom [A] 


Vdom [A 3 ] 


P [g/A 3 ] 


d [molecules/A 3 ] 


24.835 


15'318 


1.0002 


0.0334 



Table I: Length of the side of the cubic simulation domain resulting 
from equilibration, its volume, mass and molecular densities. 



much richer than what reported here. A more detailed treat- 
ment can be found in Refs. [2-5 1. In order to obtain the 
profile of the droplets, it is applied the method explained in 
Refs. EJQIlEl with AV = 1.9 x 10~ 4 [A 3 ]. Typically, the 
contact angle is extracted by superimposing a circle on the 
profile coming out from the simulations. The center and the 
radius of the circle are obtained by a fit. Here we have de- 
cided to proceed in a quite different way. Indeed, there ap- 
pears that the circle departs from the shape of the droplets in 
particular for the smallest ones where the contact angle is cal- 
culated, because its radius of curvature is weaker. We thus 
approximate the profile of the droplets by a piecewise linear 
function. The idea is to subdivide the profile into small ele- 
ments. For the points falling within an element, we calculate 
the average values for both x and y coordinates. These aver- 
age values are the usual parameters used in linear regressions. 
The contact angle is calculated from the slope of the linear 
function obtained by a linear regression on the average points 
above the contact line at most 5 [A]. The contact line point 
is assumed to be at y = 2 1//6 Oco- The contact area is simply 
given by C = nR 2 ; the contact line is L = 2%R. The overall 
interfacial surface of the dropl ets is cal culated by means of the 
formula S = 2%r 2 + 2%r(r + v V 2 — R 2 ) and their volume with 
V = (2/3)7rr 3 + 7rr 2 Vr 2 -R 2 - {n /3){r 2 - R 2 ) 3 ' 2 . 



IV. RESULTS AND DISCUSSION 
A. Water 

Table U summarizes some final results for the box of 
SPC/Fw water that is used as source in order to extract the 
droplets. Our findings are in good agreement with the original 
work for that atomistic water model [21 1. Figure [TJ shows the 
oxygen-oxygen radial distribution function. The results for 
the cumulative probability P(r) in Fig.|2]indicate that the be- 
havior of water differs slightly from that of a continuum, ho- 
mogeneous medium, except in the closest neighborhood of the 
oxygen atoms. The inset tells us that technically the first peak 
of the radial distribution function is related to the derivative 
of a step function (or Heaviside function). We thus consider a 
function of the type f(r) = (4/3)(7tr 3 /L 3 om )/(l +& i - A -''^ B ). 
Of course, this function is no more normalized to unity in the 
interval [0,Ld om ), but it provides a good approximation of the 
cumulative distribution P(r) where the first peak of the radial 
distribution function occurs. The denominator of f{r) is rem- 
iniscent of the Fermi-Dirac distribution function, which is a 
step function at low temperatures. The parameter A fixes the 
position of the first peak and the parameter B determines its 
width and height. After a simple calculation, we find for the 
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Figure 1 : Oxygen-oxygen radial distribution function of water. The 
solid line is the plot of the function «(r), Eq.[T] The parameters A and 
B are the result of a fit to the data, represented as filled circles. We 
find A = 2.7610 ± 0.0329 [A] and B = 0.0833 ± 0.0183 [A]. Inset: 
Plot of the function f(r) (see main text), with the same parameters A 
and B (solid line) compared to the same data of Fig.|2](filled circles). 
The two curves differ at most of 2.25 • 10~ 3 . 



radial distribution function the following analytic expression: 

e (A-r)/B 



u{r) = 



1 



1 

i(A-r)/B 



1 



3B l+ e ( A - r )/ B 



(1) 



This function reproduces correctly the first peak, as shown in 
Fig. [T] The main discrepancy with the data for the SPC/Fw 
model is in the neighborhood of the first minimum, corre- 
sponding to a density depletion with respect to the bulk value. 
The above method is of course unable to explain this aspect of 
the fine structure of water. 



B. Droplets 

Figure [3] shows the profile of a small droplet when approx- 
imated by a piecewise linear function. Near the contact line, 
it appears that this procedure is more precise than a circu- 
lar fit. From local averages, the resulting contact angle is al- 
ways lower (see Fig. |4j). In both cases, the macroscopic con- 
tact angle can be extracted from an average over the largest 
droplets, those of radii 60-80 [A], leading to (L = 103.1° 
(&„ = 108.1° for circular fits). Its standard deviation is 0.9% 
(0.3%) of this average value. The residual difference between 
the two methods indicates that, even for the largest droplets, 
in the close neighborhood of the contact line, the profile still 
deviates from a perfectly spherical cap. It turns out that the 
method based on circular fits underestimates the base radius. 
The average bulk molecular density is 0.0323 [molecules/A 3 ] 
(cf. Tab. [I]), with standard deviation 3.5% of it. By fitting 
the data to the modified Young equation it is extrapolated a 
value of 0oo = 111.1° (110.7° for circular fits). Especially 
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Figure 2: Cumulative probability P(r) (see main text). The dashed 
line accounts for the continuum, homogeneous behavior. We start 
from a sphere centered at the origin. When the radius of the sphere 
is larger than L ( j om /2 (see Tab. [ijl, the sphere overlaps with itself 
because of periodic boundary conditions and we can no more use 
the formula P(r) = (4/3 ^r 3 /!, 3 , . We thus employ the Monte 
Carlo principle in its elementary form. For a given radius larger than 
Lfo m /2, 30' 000 points are randomly (uniform distribution) placed in 
a cube of side L(j om and we count the fraction of them falling within 
the portions of the sphere. Inset: Magnification around the position 
of the first peak of the radial distribution function (cf. Fig.[TJ. 



for the circular profiles, we prefer the first method because 
the results indicate that around 60 [A] occurs a transition to 
the macroscopic limit (see Inset of Fig. |4] and previous ba- 
sic statistics). From a typical value for the surface tension 
of water of y = 72 [mN/m], for the line tension it is found 
K = —2.04- 10~ n [N] from circular fits and K = -5.35 • 10" 11 
[N] from local averages. This means that the forces tending 
to expand the contact area are higher according to the method 
based on local averages. In Fig. [5] we compare the contact 
area to interfacial area ratio, i.e. C/S, with the results for the 
droplets having the same final radius but with contact angle 
0oo = 108.1° (the value of 110.7° underestimates significantly 
C/S for the largest droplets). Given a final radius, the cur- 
vature of the contact line is determined numerically so as to 
have the macroscopic contact angle 9„- With respect to the 
macroscopic expectation, it is found that for the five smaller 
droplets on average the contact area expands of 4.2%, the in- 
terfacial surface shrinks of 2.8% and the volume contracts of 
7.2%. Figure [6] compares the density maps of three droplets 
of different size. These representations suggest that, for small 
droplets, a higher fraction of molecules is involved in density 
fluctuations at the interface. A simple calculation shows that 
this aspect effectively occurs if N1/N2 > pi/pi holds. The 
index 1 designates a large droplet and 2 a smaller one. N is 
the number of atoms in the fluid phase and p the bulk den- 
sity. If we compare the droplets of initial radii 30 and 60 [A], 
we find (3-15'119)/(3-l'902) > 0.0338/0.0337 = 1. Since 
even for the smallest droplet of initial radius of 20 [A] the den- 
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Figure 3: Profile of the droplet of initial radius of 30 [A]. The green 
curve represents the profile of the droplet as obtained from local av- 
erages. The straight line is a guide for the eyes and its slope amounts 
to the value used for determining the contact angle. Inset: Profile 
from a circular fit to the data in the neighborhood of contact line. 
The straight line is the tangent at the contact-line point. 



sity of water is close to the bulk value and the surface thick- 
ness is around 3 [A], we conclude that for sure larger droplets 
have a reduced fraction of molecules at the interface. In other 
words, as the droplet increases in size, its interface grows and 
the fraction of molecules fluctuating at the interface becomes 
smaller. The droplet of initial radius 60 [A] was also simulated 
at different temperatures up to 450 [K]. For this temperature 
increase, it is found that the contact angle varies with good 
approximation linearly, as well as the molecular density. In 
contrast to what reported in Ref. [25 1 the contact angle varies 
over this temperature range only of a few degrees. Preliminary 
results using coarse-grained models confirm our trend. 



V. CONCLUSIONS 

Contact angle measurements and comparisons with the pre- 
dictions of Young equation are generally carried out under the 
assumption of a spherical shape of droplets. Figure [5] shows 
that, when the contribution of the corrective term to Young 
equation and/or other size effects J6] [TTl |24] is more impor- 
tant, the droplets deviate significantly from the macroscopic 
expectation. Furthermore, below the initial radius of 30 [A], 
the contact angle can no more be derived accurately from the 
tangent to a circular profile at the contact line (see Figs.[3]and 
0). Actually, from our analysis it clearly emerges that, for 
small droplets, fluctuations of the surface thickness are of ma- 
jor relevance, resulting in the deformation of their spherical 
shape. For small droplets the effect is more marked presum- 
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Figure 4: Contact angle dependence on the final radius from the 
two approximations of the droplet profiles (squares for circular fits 
and circles for local averages). Inset: Fitting of the data to modified 
Young equation with contact angles resulting from circular fits, Top, 
and from local averages, Bottom. 
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Figure 5: Contact area to overall interfacial surface ratio as a func- 
tion of the final radius of the droplets. The straight line is the macro- 
scopic expectation under the assumption of a perfectly spherical pro- 
file. 



ably because of their reduced size and the short-range nature 
of non-bonded interactions (cf. Ref. [24]): cohesive forces 
near the contact line are weaker and the contact area would 
tend to expand, leading to lower contact angles. Interestingly, 
it also appears that the macroscopic limit is not reached grad- 
ually, but with a possible transition around the initial radius 
of 60 [A]. On the other hand, the predictions of the modi- 
fied Young equation are more precisely recovered for initial 
radii around 40 [A]. Finally, the SPC/Fw model for water is 
extensively validated j2Tl and the Lennard-Jones parameter 
£co used here resulted from a previous calibration accord- 
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Figure 6: Density maps for the droplets of initial radii 30, 60 and 80 [A], in this order. Color code based on molecular density. The thickness 
of the interface is in all cases around 3 [A]. Concerning the first map, the horizontal shift relative to the representation of Figgis due to the 
fact that the first radial bin coincides with zero and the coarser binning necessary for this kind of representation. 



ing to recent measurements carried out under the most ideal 
conditions E3l [24). Our macroscopic contact angle 0oo of 
111.1° (circular fits and extrapolated value) underestimates 
the experimental value of 127° 11231 . taken as reference for 
calibration |24|. We ascribe this discrepancy mostly to the 
different cutoff scheme and the inclusion of long-range inter- 
actions 11261 . Clearly, our analysis of profiles and the main 
conclusions drawn for the smallest droplets in the hydropho- 
bic regime are robust under small changes of the simulation 
settings. Indeed, the profiles of these droplets will be affected 
to a larger extent by the fluctuations of the surface thickness 
and in turn no more fully spherical. The accuracy implied by 
the various simulation settings confer a substantial degree of 
confidence to the findings reported here. 



VI. NOMENCLATURE 



A,B 


A 


model parameters 


C 


A 2 


contact area of droplets 


d 


A- 3 


molecular density 


f,u 




model functions 


S 




radial distribution function 


L 


A 


length of contact line of droplets 


^dom 


A 


side length of simulation domain 


/ 


A 


carbon bond length in graphene 


m 


g-mol -1 


mass of atoms 


N 




number of atoms 


n 


A- 1 


number of spherical shells per 


length 


P 




cumulative radial distribution for 






gen atoms 



R 


A 


r 


A 


Ar 


A 


S 


A 2 


V 


A 3 


Vdom 


A 3 


w 




x,y,h 


A 


e 


kcalmol 


7 


N-nr 1 


K 


N 


<3 


A 


P 


g-cm -3 


e 





-1 



base radius of droplets 

radius of droplets 

width of a spherical shell 

overall interfacial surface of droplets 

volume of droplets 

volume of simulation domain 

average number of oxygen atoms in a 

spherical shell 

cartesian coordinates 

interaction parameter for Lennard-Jones 

potential 

surface/interfacial tension 
line tension 

interaction parameter for Lennard-Jones 
potential 
mass density 
contact angle 
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